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ABSTRACT 

We utilize a chemical equilibrium code in order to study the condensation process 
which occurs in protoplanetary discs during the formation of the first solids. The model 
specifically focuses on the thermodynamic behaviour on the solid species assuming the 
regular solution model. For each solution, we establish the relationship between the 
activity of the species, the composition and the temperature using experimental data 
from the literature. We then apply the Gibbs free energy minimization method and 
study the resulting condensation sequence for a range of temperatures and pressures 
within a protoplanetary disc. 

Our results using the regular solution model show that grains condense over a large 
temperature range and therefore throughout a large portion of the disc. In the high 
temperature region (T ^ 1400 K) hibonite and gehlenite dominate and we find that 
the formation of corundum is sensitive to the pressure. The mid-temperature region 
is dominated by Fe( s ) and silicates such as Mg 2 Si0 4 and MgSi0 3 . The chemistry 
of forsterite and enstatite are strictly related, and our simulations show a sequence 
of forsterite-enstatite-forsterite with decreasing temperature and the abundance of 
the first high temperature peak of forsterite is also pressure sensitive. In the low 
temperature regions (T ^ 600 K) a range of iron compounds (FeS, Fe 2 Si0 3 , FeAl 2 3 ) 
form. We find that all the condensation sequences move towards lower temperature as 
the pressure decreases. 

We also run simulations using the ideal solution model and see clear differences in 
the resulting condensation sequences with changing solution model. In particular, we 
find that the turning point in which forsterite replaces enstatite in the low temperature 
region is sensitive to the solution model. In this same temperature region, fayalite is 
the most stable compound for the regular solution, while magnetite replaces fayalite 
in the ideal solution model at the lowest values of temperature. 

Our results show that the ideal solution model is often a poor approximation to 
experimental data at most temperatures important in protoplanetary discs. We find 
some important differences in the resulting condensation sequences when using the 
regular solution model, and suggest that this model should provide a more realistic 
condensation sequence. 

Key words: astrochemistry, planetary systems: protoplanetary discs, star: pre-main 
sequence, circumstellar matter 



1 INTRODUCTION 



Studying protoplanetary discs, and more particularly their 
dust content, is crucial for understanding the creation of 
* E-mail: fpignatale@astro.swin.edu.au the solar system. Solar system bodies such as planets, me- 
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teorites and comets are all created from solids during the 
protoplanetary disc phase. 

The quantity and the quality of observations of proto- 
planetary discs around young stars has increased substan- 
tially in the last years: Spitzer Space Telescope has provided 
a large amount of useful data whose interpretation provides 
information on the composition, growth and thermal pro- 
cessing of dust grains in discs dOlofsson et al.||2009| iMeeusI 
et al.|2009||Bouwman et al.|2008||Kessler-Silacci et al.|2006[ 



Geers et ai.||2006[ ). The processes of grain formation are in- 
trinsically related to the chemistry of gas and dust. Indeed, 
the material in the discs is a mixture of condensates with 



different chemical composition and properties (Henning & 
Meeus|2009l |Bergin|2009l >. 



The study of objects in our solar system, including as- 
teroids, meteorites and comets provides more direct and ac- 
cessible evidence of the chemistry of the young solar system. 
Meteorites are characterized by heterogeneous compositions 



( Krot et al. 2005 1 and petrographic and mineral analysis 



show mixtures of compounds that are the result of conden- 
sation processes in different environments and times during 



the protoplanetary disc phase ( |Scott fc Krot 20051. Even 
comets, among the oldest and most distant objects of our 
solar system, show this heterogeneity ( Brownlee|2005 1 . The 
Stardust mission (Brownlee et al. 2003| ) returned samples 



from comet Wild 2 and studies of their composition ( Naka- 
mura et al. 2008 Zolcns ky et al.|[2006[ ) show complex ther- 



modynamic and dynamic scenarios where particles proba- 



bly formed in the inner regions of the solar system ( Liffman 



20061 which experienced high temperature processes were 



found close to other material produced by cold thermody- 
namics. 

This paper is the first in a series in which we aim to 
determine the composition of both the gas and solid phases 
in various regions of protoplanetary discs in order to un- 
derstand their chemistry and study the dust content. We 
assume thermodynamic equilibrium and solve for the con- 



densation sequence using the regular solution model ( DeHoff 
|1993[ ) for the behaviour of solid solutions. In this paper we 
present the regular solution model and compare the conden- 
sation results with those obtained using an ideal solution for 
a range of temperatures at fixed pressures. In our next pa- 
pers, we will apply this technique to specific ID and 2D discs 
models. 

The Gibbs energy minimization technique for equilib- 
rium calculation is a powerful method for understanding the 
evolution of chemistry in a complex system ( |DeHoff||1993[ ). 
However, the predictions using this technique are sensitive 
to providing rigorous mathematical descriptions of solution 
behaviour, in particular, providing allowances for the effect 
of concentration and temperature on the activity of species 
dissolved into a solid state solution (DeHoff 1993). In this 



study, we include the effect of temperature of solution be- 
haviour by using the regular solution model, which is su- 
perior to the assumption used in previous studies in which 
ideal mixing is assumed for many of the solid phases and 
the effect of temperature on solution behaviour is ignored 
(e.g. |Pasek et al.|2005||Gail|1998[ ). It is well established that 
these effects are important in complex phase equilibria. 

The outline of this paper is as follows: in section 2 we 
describe the physical and chemical conditions assumed to be 
present in protoplanetary discs. We also introduce the Gibbs 



free energy minimization method to solve for the condensa- 
tion sequence and the steps followed to determine our final 
thermodynamic model. In section 3 we discuss the behavior 
of solid solutions and present a range of activity-composition 
relations derived from laboratory experiments. In section 4, 
we present the results of our simulations and we compare 
them with the ideal solution models (section 5) and with 
previous modeling work and recent observations (section 6). 
Conclusions are presented in section 7. 



2 THE MODEL 

Protoplanetary discs results from conservation of the angu- 
lar momentum of a collapsing, rotating cloud of gas and 
dust. Detailed studies have been made to define the struc- 
ture of protoplanetary di scs (|Brauer et aT7||2008| |D'Alessio| 



et al. 2004 1999 1998 Gail 2004 2001 1998) and the 



physics involved in the evolution of the resulting gas-dust 
system ( jBirnstiel et al.|2010|[bullemond et al.|2008[ ). In this 
work, we calculate the condensation sequence for an initial 
chemistry over a range of temperatures and at fixed values 
of pressure which are appropriate for protoplanetary discs. 

In this section, we describe the conditions within our 
model, the chemical equilibrium code used to study the con- 
densation sequence, along with assumptions and limitations 
of the code. 



2.1 Physical properties of the disc 

Different disc models return different temperature and pres- 
sure distributions. These can be described as functions of 
the disc radius (r) in ID models and both the radius and 
the distance from the mid-plane (r,z) for 2D models. In this 
model, we do not subscribe to any specific disc model, but 
instead chosen range of temperatures and specific pressure 
values that are relevant in protoplanetary discs. These val- 
ues are not meant to represent specific locations in the disc 
but to constrain the values of temperature and pressure we 
use in the our simulations and define three main regions 
(high, middle, low temperatures) in the disc in which con- 
densation occurs. Furthermore, condensation sequences cal- 
culated over a range of temperatures at fixed values of pres- 
sure allow us to compare our results with previous work and 
thermodynamical predictions. 

For this work, we use the disc model of D'Alcssio et al. 



11998) to guide us to find suitable values of T and P. 



D'Alessio et al. (19981 studied the vertical structure, tem- 



perature, and surface density of a disc whose heating pro- 
cesses include viscous dissipation, radioactive decay, cosmic 
rays, and stellar irradiation. They chose a typical T Tauri 
star with M* = 0.5 Mq, R* = 2R , T* = 4000 K, and 
M = 10 _8 M s yr _1 . 

We assume an irradiated disc with a constant accre- 
tion rate and consider only the inner 5 AU. Most of the 
dust mass is located near the mid-plane of the disc, so the 
study of temperature and pressure was done at the mid- 
plane. |D'Alessio et al.| ( |1998[ ) calculated the radial temper- 
ature profile at the mid-plane, the surface density, and the 
different vertical heights above the mid-plane. The tempera- 
ture at the mid-plane, T c , ranges from 60 to 1900 K for radii 
between 5 and 0.05 AU. The values of P and T c are shown in 
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Table 1. Temperature and pressure values for three typical radii 
calculated from the model oflD'Alessio et al.l (|1998|l. 



Parameters 


r = 0.05 AU 


r = 1 


AU 


r = 5AU 


Tc (K) 


1900 


240 




60 


P (bars) 


2.1 x 10~ 3 


2.2 X 


io- 7 


9.4 x 10" 9 



Table [T] for three typical disc radii: r = 0.05, 1.0 and 5 AU. 
In the condensation studies we present in this work, we use 
a temperature range from 50 to 1850 K for fixed pressures 



10 _J ,10" 



10 -8 bar. 



2.2 Gibbs free energy minimisation 

Our condensation model of protoplanetary discs is based on 
the determination of the equilibrium of the initial system 
given a set of temperatures, pressures and initial composi- 
tion. For transformations occurring at constant temperature 
and pressure, the relative stability of the system is deter- 
mined by its Gibbs free energy. The general definition of 
Gibbs free energy is 



G = H — TS , 



(1) 



where H is the enthalpy, T the absolute temperature and S 
the entropy. The laws of classical thermodynamics determine 
that the system will be in its most stable state if the Gibbs 
free energy is minimized (DeHoff 1993; Porter & Easterling 
19921. In a system with multiple species, the total Gibbs 



free energy, Gt, is given by 



(2) 



where Xi is the number of moles of species i, and Gi is the 
Gibbs free energy of species i, and N is the total number of 
possible species defined for the system. Gi is a function of 
the temperature, the mole fraction and solution parameters 



G l = G° i + RT\n<y t + RTlnXi 



(3) 



where G; is the standard Gibbs free energy of formation of 
the species i, R is the gas constant, T the temperature, Xi is 
the molar fraction of the species in the solution phase, and 7, 
is the activity coefficient of species i. The variation of fi with 
temperature and composition is determined experimentally. 

The determination of the equilibrium is equivalent to 
finding the set of values which minimizes the function in 
equation |2| and satisfies the mass balance constraint given 
by 



= bj {j = 1,2, ...,m) . 



(4) 



where there are in different types of elements and CLij is the 
number of atoms of element j in species i, Xi is the number of 
moles of species i (such that the molar fraction Xi = Xi/xt, 
where xt is the total number of moles), and bj is the total 
number of moles of element j. 

Therefore, the parameters required to solve equation [3] 
include the temperature, T, the number of moles Xi of each 
of the TV species i and the total number of moles x t in the 
initial composition of the system, the standard Gibbs free 



energy of formation of each species i, Gi, and the activity 
coefficients ji, which depend on the behaviour of the species 
i in the solution, for a given temperature, T. 

The Gibbs energy minimisation method is a technique 
widely used for geology, chemistry, astrophysics, metallurgy, 
materials and chemical engineering and industrial purposes 
for understa nding comp lex equilibrium calculation at high 
temperature (Hack 1996 1. For our purpose, we use the HSCQ 
software package (Roine 20071, which uses the Gibbs free 



energy minimization method of White et al. ( 19581. 

The thermodynamic data for each compound in our list 
are taken from the database provided by HSC. The list of 
references is very long and not reported here. For each com- 
pound the database reports the enthalpy, Ft, entropy, S, and 
heat capacity, C, along with a range of other important pa- 



rameters. We refer the reader to Roine (20071 and references 



within, plus the HSC Database module for full details. 



2.3 Method and limitations 

In the previous section, we introduced the concept of equilib- 
rium which assumes that the system has reached the most 
stable state for a given T, P and composition. This is not 
completely true for a protoplanetary disc. For example, the 
existence of non-crystalline material in protoplanetary discs 
is a sign that the disc is not in equilibrium and at low tem- 
peratures reaction rates and transport phenomena can in- 
deed be so slow as to make equilibrium unachievable. How- 
ever it is assumed that the lifetime of a disc (~6-30 Myrs) 



(Carpenter et al. 2005 Haisch et al. 20011 is long enough 



to overcome the kinetic barriers such as mass transfer of 
species, reaction rates and surface phenomena, thus allow- 
ing the disc to reach an equilibrium state. To ppani et al.| 
( |2006[ ) studied the reaction rates of many compounds in our 
list (e.g. forsterite, anorthite, spinel, corundum, pyroxene, 
melilite) also focusing on the influence of time, tempera- 
ture and gas composition in the reaction rates. Their results 
show that a mean time of 1 hour is required to return sta- 
ble phases. Furthermore, computational techniques based on 
the assumption of equilibrium have been successfully used 
to study and understand reactions in metallurgical systems 
and in the Earth's geology processes ( Belov|2 002, and refer- 
ences therein). 

We assume that the system is only composed of solid 
and gas states. While the condensation model developed 



by Yoneda & Grossman ( 1995 \ showed that liquids can be 



formed at atmospheric pressures, the typical pressure of a 
protoplanetary disc is much lower than this (P ^ 10 -2 bar 
- see Table [TJ . The absence of liquid phase in this pressure 
regime can also be seen in the T-P phase diagrams of com- 
pounds in our list, and hence we ignore the liquid phase. 

The range of possible species produced in the model are 
derived from a limited number of elements. We choose the 
fifteen most abundant elements of the Sun taken from Pasck 



et al. (20051, who use the solar photosphere abundance ob- 



servations of Grevesse & Sauval ( 1998 1, Allende Prieto et al 
p002] > and I Allende Prieto et aL] ( |2001| ) (see Table |2j) from a 



total amount equal to 100 kmol. We assume that the Sun's 



HSC website: http://www.outotec.com 



4 Pignatale et al. 



photosphere has changed little over the time and so the cur- 
rent composition can safely be used for the initial conditions 
in our model. 

The chondrites, which formed early in the history of the 
solar system, have the same abundance for most elements 
as the photosphere (Scott & Krot 20061, which supports 



Table 2. Abundance of elements in our system, taken from |Pasek| 
et al. (2005), for a total amount equal to 100 kmol. 



our assumption of a relatively constant stellar composition. 



(While Asplund et al. ( 2009 \ have a revised photospheric 



composition, we use the values of Pasek et al. ( 2005 \ to 



compare our results to previous work.) 

The potential number of species that can be formed 
from the combination of these 15 elements considered in 
this system is over 700. Our aim is to improve upon previ- 
ous works in this field by using the regular solution model, 
and for the non-ideal case the HSC software is limited in the 
number of species it can handle. Using the 15 elements in 
Table[2] we follow [Pase"k et al.|prjQ5| to prepare an initial list 
of species likely to form. Using the tables in Appendix Al, 
we first remove all compounds that include the extra 8 solar 
elements which we exclude (i.e. Cr, P, Mn, CI, K, Ti, Co, F). 
We next use this list and assume ideal solution behaviour 
for just two phases, gas and solids, which are assumed to be 
pure phases. We run this ideal model for temperatures rang- 
ing from 50-1850 K for fixed pressures 10~ 3 , 10~ 4 , 10~ 8 
bar. Species which did not appear in the system within the 
T and P ranges of our models were deleted, such as carbides, 
nitrides and some sulfide species. This procedure eliminated 
species that are very unlikely to be formed, while allow- 
ing the software to run efficiently. However, we kept some 
species, such as solid Al and various oxides which are clearly 
important in the chemistry of the condensation sequence of 
the solar system and other protoplanetary discs, to verify if 
these compounds, for our range of temperatures and pres- 
sures, are involved in the formation of more complex species 
or not. 

Carbon is only present in our system in the main 
gaseous compounds together with solid graphite. Ideal cal- 
culations made during the definition of our initial system 
show that carbon solid compounds are very unlikely to form 
for our range of pressures and temperatures in equilibrium 
conditions. The majority of the initial C(g) is involved in 
the formation of gas compounds like CO and CH 4 (see sec- 
tion j4|, due to the presence of large amounts of oxygen and 
hydrogen in our system. Indeed the Gibbs free energies of 
these gas compounds makes carbon gases the most stable 
carbon compounds. Furthermore, the formation of graphite 
is strictly related to the environmental condition in which a 



large amount of long chains of carbon gas is found (Pasek 
et al.|[2005| . Our equilibrium calculations show that lower- 
ing the temperature results in C(g) being replaced by CO(g) 
and there is not enough C(g) in the range of condensation 
temperature available to form graphite. 

The complete list of species and phases used in our 
model is shown in Table [3] 



3 BEHAVIOUR OF SOLID SOLUTIONS 

We assume that solid species can mix together to form solid 
solutions, as these phases are observed on Earth as miner- 
als and in meteorites. The choice of solid solutions is jus- 
tified from previous thermodynamic studies. The behaviour 



Element 


Abundance (kmol) 


ii 


Q 1 


He 


8.89 





4.46 x lO" 2 


c 


2.23 x 10~ 2 


Ne 


1.09 x 10~ 2 


N 


7.57 x 10~ 3 


Mg 


3.46 x 10~ 3 


Si 


3.30 x 10~ 3 


Fe 


2.88 x IO- 3 


S 


1.44 x 10~ 3 


Al 


2.81 x 10~ 4 


Ar 


2.29 x 10~ 4 


Ca 


2.04 x 10~ 4 


Ka 


1.90 x 10~ 4 


Ni 


1.62 x 10" 4 



of each solid solution is studied, in order to determine their 
relation between composition, temperature and activity pa- 
rameter. 

The choice of the solution behaviour which we use in 
this study is based on both laboratory and theoretical ther- 
modynamic studies for each phase found in the literature. 
Most of these works focused on the relations between species 
composition and activity for one temperature, whereas ob- 
servations of activities clearly show non-ideal solution be- 
haviour evolving with temperature. The solution behaviours 
can be quite complex with several parameters involved in 
the reactions. There is no comprehensive scientific treatment 
of solution behaviour and no single approach to modelling 
multicomponent solutions, but there are several developed 
models available (DcHoff 19931. We choose the regular so- 
lution model to establish relations between the temperature 
in which the experiments were made, called reference tem- 
perature, the composition of the solution, a range of tem- 
peratures, the activity and the activity coefficient describing 
a particular compound in its phase. 

Whilst the thermodynamics of pure species is largely 
uncontroversial, solution behaviour is difficult to mea- 
sure and model, particularly for complex oxide and sul- 
fide phases. Therefore critically examining the solution be- 
haviour of species important for understanding protoplan- 
etary dust formation is vital in developing a rigorous un- 
derstanding of the system. In particular, we have sourced 
the best available solution data from geological and metal- 
lurgical studies to incorporate into studies of astrophysical 
dust, as the underlying physical chemistry of these different 
systems is fundamentally the same. In this work we have at- 
tempted to include the best available solution behaviour and 
where we have made assumptions, we state these clearly as 
a means by which to identify experimental conditions that 
need to be studied in the laboratory and stimulate the de- 
velopment of approximate solution models. 



3.1 The regular solution model 

The regular solution model is based on two fundamental 
assumptions ( DeHoff|1993 1: (i) the entropy of mixing is the 
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Table 3. List of phases with their species used in our model. 



Gas phase Solid phases 



Al 


Mg 


SiS Olivine 


Melilite 


Magnesiowustite 


Pure phases 


C 


N 


Mg 2 Si0 4 


Ca 2 Al 2 Si0 7 


FeO 


Na 2 Si0 3 


CH 4 


N 2 


Fe 2 Si0 4 


Ca 2 MgSi 2 7 


MgO 


CaAl 12 O lg 


CO 


NH 3 








H a O 


Ca 


Na 


Orthopyroxene 


Fassaite 


Sulfide 


A1 2 3 


CaO 


Ne 


MgSiO a 


CaAl 2 Si0 6 


FeS 


CaO 


Fe 


Ni 


FeSi0 3 


CaMgSi 2 6 


NiS 


Fe 2 3 


FeO 


O 




Ni 3 S 2 


Fe 3 4 


H 


o 2 








Si0 2 


H 2 


OH 


Plagioclase 


Spinel 


Metal 


C 


H a O 


s 


CaAl 2 Si 2 O g 


MgAl 2 4 


Fe 




HS 


Si 


NaAlSi 3 Og 


FeAl 2 4 


Ni 




H 2 S 


SiO 






Si 




He 


SiQ 2 






Al 





same as that for an ideal solution, and (ii) the enthalpy of 
solution is not zero, as in ideal solution, but is a function 
of composition. As such, the model gives a simple relation 
between the activity coefficient and the temperature: if 7^ 
is the activity coefficient of one species at a temperature Ti 
and 7t 2 is the activity coefficient at another temperature T 2 
(for the same composition) then 



ln_7Tj_ 
ln7 T2 



T2 
2i 



7Ti 



T2/T1 
It* > 



(5) 



(6) 



In the absence of activity data for a species (k) , a simple re- 
lation between the temperature and the activity coefficients 
was introduced, based on the observations of the evolution 
of the behaviour of the species with temperature given by 

0.99 To/T , 



7* cn 



(7) 



where To is the reference temperature for the study of the 
phase. Thus, the activity coefficient is near the ideal (i.e. 
7 ~ 1) for high temperatures when T ^ To, while it moves 
away from the ideal as the temperature drops (i.e. T ^ To). 
This assumption provides a reasonable method to extrapo- 
late existing solution data to conditions where there are no 
measurements of solution behaviour, and is tested in sec- 
tion [3]3] 

Most thermodynamic literature report experimental 
values of activity at several values of composition and To. 
The activity, a, and the activity coefficient, 7, are linked by 
the relation 



a>k(Xk, T ) = 7feX fe , 



(8) 



where X is the composition of the studied species fc. The 
relations in equations |8| and |6| allow us to derive a poly- 
nomial for the activity coefficient given by 

Tk(X ki T) = a hi^M = 7fc(Xfc)To) (Tom (9) 

The meaning of activity and activity coefficient can be 
understood from equation |S|. If 7 = 1 (ideal) the activity 
of component k is equal to its mole fraction Xj. and the be- 
havior of k is completely determined by its composition. If 



the value of 7^ 7^ 1 (as in the regular solution model) the 
component acts as if the solution contains more (7^ > 1) 
or less (jk < 1) than the mole fraction X suggests. Further- 
more, from equation Q the activity coefficient 7 provides 
information about the departure from the ideal behavior. 

The regular solution model is a quantitative explana- 
tion of the non-ideal behaviour of solid solution, but there 
are some limitations. In a dilute solution of component b 
(solute) in the component a (solvent), Xb approaches or 
X a approaches 1. In this limit the activity coefficient for the 
solute component 76 becomes close to the Henry's law co- 
efficient. A similar behaviour for this kind of solution can 
be seen in the opposite limit of the composition range when 
component a is dilute in component b. This symmetry in 
the extreme values of composition implies that the proper- 
ties associated with an atom of component b surrounded by 
atoms of component a are the same as the opposite case in 
which an atom of component a is dilute in a component 6 
HDeHoff|1993] ). Clearly this is not the case. Furthermore, the 
regular solution model returns a parabolic AH mlx (heat of 
mixing) function that is also symmetrical with respect to 
composition, whereas it is often unsymmetrical, especially 



in liquid solutions ( Porter fe Easterling|1992| ). 



3.2 Solid solutions and their behaviour 

For each solution and each species, using thermodynamic 
laboratory experiments and theoretical models from the lit- 
erature, and the use of the regular solution model as de- 
scribed above, we can now establish a relationship between 
the composition, temperature and activity coefficient of the 
nine solid phases listed in Table [3] 

As a general principle when choosing the data for the 
activity, we have tried to find direct experimental data re- 
porting activity values of compounds. These data are pre- 
ferred over indirect methods derived via mathematical tech- 
niques, such as the Margules and Gibbs Duhem equations, 
since these techniques inherently approximate solutions (see 
e.g. |Bodsworth fe Appleton||I5o5l |Anderson|[2005l |Darken| 
& Gurry 1953 1 which could affect the results. Thus where 



possible we preferred to choose data derived from direct lab- 
oratory experiments. 
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FeO Composition 




0.2 0.4 0.6 0.8 
NaAISi 3 O s Composition 




Fe Composition Ni Composition 

Figure 1. Activity-composition functions for various phases. The 
data points are from experimental data and the curves are our 
polynomial fit to the data. The straight line is the ideal solution. 
Top (a): Magnesio-wustite solution: activity-composition func- 
tions for_FeO (left) and MgO (right) at T = 1573.15 K. Data 
from|Nafzigcr (1973). Middle (b): Plagioclase solution: activity- 
composition functions for albite ( NaAlSi 3 O s — le ft) and anorthitc 
(CaAl 2 Si 2 8 - right). Data from |Newton et al.| ( |1980| . Bottom 
(c): Metal solution: activity-composition functions for Fe and Ni 
at T = 1273.15 K. Data from |Conard et a.1 ■ | ( fl979| l ■ 



3.2.2 Sulfide phase 

The sulfide phase is composed of FeS (troilite), NiS, and 
Ni 3 S 2 . The activity data of FeS and Ni 3 S 2 was studied by 



Byerley & Takebe (19721, who measured the sulfur activ- 



ity for the Fe-Ni-S solid ternary system in equilibrium with 
gaseous mixtures of H 2 S and H 2 at To = 1523.15 K, and us- 
ing the Gibbs-Duhem equation. The activity data for these 
two species were presented in an isoactivity Fe-Ni-S ternary 
diagram in Byerley & Takebe ( 1972| ) . An analytical function 
of composition for the activity coefficient could not be de- 
duced, and hence the determination of the activity had to 
be made using the composition values of sulfide compounds 
found in our first test simulations using the ideal solution 
model. 

For most temperatures in our calculations where FeS 
and Ni 3 S 2 appear, TVfcs = 1.43xl0" 3 kmol, iV N i 3 s 2 = 6.13x 
10" 5 kmol, and N ms = 1.5 x 10" 6 kmol, where N is the 
number of mols. Thus the composition of Fe, S and Ni in 
the system are X s = 0.49, X Fe = 0.46 and X Ni — 0.05. 

Plotting these data on the isoactivity diagram of |By-| 



erley & Takebe (19721, the activity was calculated to be 

0.93, with a composition rate 
= 0.04, and 7 Fe s = 0.97 and 
Thus 



dNi 3 s 2 = 0.01 and a FeS = 
of X FeS = 0.96 and -X"Ni 3 s^ 
7ni,s, = 0.25 at T = 1523.15 K 



7FeS(T) = 

7Ni.s. (T) = 



0.97 
0.25 



T /T 
T /T 



(12) 
(13) 



As no activity data could be found in literature for NiS, 
equation |7f provides the following relationship 

lms(T) =0.99 To/T , (14) 

where T = 1523.15 K. 



3.2.1 Magnesio-wustite phase 

MgO and FeO can mix together to form the magnesio- 
wustite phase. Nafziger ( 1973 I studied the behaviour of 
these compounds in the MgO-FeO-F 2 O s system at To = 
1573.15 K, providing activities of both monoxides according 
their composition in the system. 

Using these experimental data, an explicit relation be- 
tween the activity and the composition for FeO and MgO 
can be deduced. We fit a polynomial function to the data 
using a least squares fit algorithm with weight and standard 
deviation to determine the best polynomial order. Thus an 
activity coefficient-composition relation could be derived for 
To = 1573.15 K and then for the entire temperature range 
using the regular solution model (equation |9| given by 



lFe o(T,X Fe0 ) = (-10.62Xf. eO + 38.59A:| e0 - 54.48X| eD 



T /T 



+38.86X^0 - 15.43X Fe0 + 4.07) 

r yMgo(T,XM g o) = (— 23.36A^ /s0 + 77.79Xm 9 o 
-102.20Al fgO + 67.88X 2 Mg0 - 24.20X MgO 
+5.09) To/T . 

These relations are shown in Fig. [la. 



(10) 



(11) 



3.2.3 Olivine phase 

This phase is composed of forsterite (Mg 2 Si0 4 ) and fayalite 
(Fe 2 Si0 4 ). The study of the activity of this phase is based 



on data from Nafziger & Muan (19671, who experimentally 



determined the activity of olivine species at To = 1473.15 K 
for the MgO-FeO-Si0 2 system. 

Depending on the initial composition of MgO, FeO, and 
Si0 2 , five solid solutions can exist - see Figure 1 of .Nafziger 
& Muan ( 19671. Assuming ideal solution behaviour, our cal- 



culations show that olivine and pyroxene species coexist for 
the range of temperature in our system. So in this study we 
chose to focus on the olivine + pyroxene case. 

Pyroxene and olivine species are produced by the fol- 
lowing reactions 

Fe + Si0 2 + 1 2 — > FeSi0 3 (pyroxenephase), 
2 Fe + Si0 2 + 2 — > Fe 2 Si0 4 (olivinephase), 

and are in equilibrium following the reaction: 

2FeSi0 3 + Mg 2 Si0 4 ^ Fe 2 Si0 4 + 2MgSi0 3 . 

Setting the oxygen partial pressure and enstatite com- 
position, and using relations between the standard free en- 



ergy of the reaction, activity, and composition (Nafziger 
& Muan 1967), the activity-composition of olivine species 
could be deduced. 
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With these data we are able to extract activities for 
Fe 2 Si0 4 and Mg 2 Si0 4 using the relations 



for albite (Ab) and anorthite (An) respectively. Newton 



aFc 2 Si0 4 = a FcS i . 5 o 2 > 



ClMg 2 Si0 4 



a MgSi, 



(15) 
(16) 



presented by Anderson & Crerar (19931. Fitting a poly- 



nomial to the experimental data presented by |Nafziger fc| 
|Muan| ( |1967[ ) for the olivine phase, an activity coefficient- 
composition relation for the two species could be deduced, 
firstly for To = 1473.15 K and then for all temperatures, 
given by 



"fFay(T, X Fay ) = 
"lFo(T, X Fo ) = 



(3.80X Fay 
-0A5X Fay + 0.24) To/T 
(7.58X Fo 
-5.19X Fo + 0.88) To/T 



A9X Fay + 5.6QX 2 Fay 



(17) 



19.19X Fo + 16.91X Fo 



(18) 



where Fay is the fayalite (Fe 2 Si0 4 ) and Fo the forsterite 
(Mg 2 Si0 4 ). 



3.2.4 Orthopyroxene phase 

The activities of enstatite (MgSi0 3 ) and ferrosilite (FeSi0 3 ) 
in orthopyroxene solution were determined experimentally 
by Saxena & Ghose (19711. They assume that the silicate 



framework does not change significantly and its contribution 
to the change in activity with temperature may be neglected, 
and therefore ctp° px ' is similar to ap^sioj ■ 

The most complete data from their study was for To = 
873.15 K, and we used these data to establish our activity 
coefficient as a function of temperature and composition. 

From the experimental data curves, we fitted a poly- 
nomial and using the regular solution model, the activity 
coefficients as a function of temperature and composition 
for the enstatite (En) and ferrosilite (Fr) were derived as 
follows: 



7En(T,X En ) = (12.90X|„ - 35.62X| n + 35.6l| n 

-15.54 Fn + 3.63) T ° /T , (19) 

7Fr(T, X Fr ) = (5.72X Fr - 13.04X Fr + 10.44X Fr 

-4.12X Fr + 2.01) T ° /T , (20) 

where T = 873.15 K. 



3.2.5 Plagioclase phase 

The activities of anorthite (CaAl 2 Si 2 O g ) and albite 
(NaAlSi 3 O g ) of the plagioclase solution must be treated 



with a non- ideal model, as shown by Newton et al. (19801. 
|Kerrick fc Darken| ( |1975[ ) derived an ideal model for pla- 
gioclase solution from statistical mechanical considerations, 
based on the aluminium avoidance principle (Al-O-Al bonds 
in CaAl 2 Si 2 O g and NaAlSi 3 O g species are believed energet- 
ically unfavorable). They deduced activity-composition re- 
lationships 



aAb(X) = X 2 Ab (2 — X Ab) , 

aAn(X) = -Xau(1- + Xau) 2 , 



(21) 
(22) 



et al. ( 1980 I used experimental studies of plagioclase to 



show that the solid solution of Kerrick & Darken ( 1975 1 



does not follow ideal behaviour. Specifically, they found that 
the measured enthalpy of solution is non-zero, contrary to 



ideal model assumption. Newton et al. ( 1980 \ found activ- 



ity coefficient-composition relationships from the measured 
excess enthalpy of mixing to be given by 

j Ab = exp[X An (W A b + 2(WAn-WAb)XAb)/RT], (23) 
7 a„ = exp[Xl b (WAn + 2(W Ab -W An )XAn.)/RT], (24) 

where W An and W A b are the Margules parameters (which 
are constants at constant temperature and pressure), and 
X An and X A b are respectively the mol fractions of anor- 



thite and albite. In their system, Newton et al. ( 1980) found 
W An = 2.025 Kcal and W Ab = 6.746 Kcal. 

In order to derive the activity and then the activity 
coefficient polynomials using equation jjy, these coefficient 
must be multiplied by the composition functions for activity 
found by |Kerrick fc Darken] ( |1975[ ) (equations (2l| and p2|), 
and thus the theoretical activity coefficients (from 7 = a/X) 
which are required by the HSC software were found to be: 

-f Ab (T, X Ab ) = (O.S7X 3 Ab - 1.57X 2 Ab + lA5X Ab 

+0.24) To/T , (25) 
1An (T,X An ) = (-0.55Xl n + 1.23X A n+0.33) To/T ,(26) 

where To = 1000 K. These relations are shown in Fig. [TJs. 



3.2.6 Melilite phase 

In order to determine the behaviour of the melilite phase 
and derive the activity-composition relations of gehlenite 
(Ca 2 Al 2 Si0 7 ) and akermanite (Ca 2 MgSi 2 7 ) in solid solu- 
tion, we used the work of Charlu et al. ( 1981 1 and Walbaum 
( |1973[ ) and the analysis of the melilite system from Beckett 
qi994p . We also made use of the MELTS Supplemental Cal- 



culator described by Ghiorso & Sack ( 1995 1 and Asimow & 
|Ghiorso| ( |1998 1 and references therein in order to derive the 
activity-composition data for gehlenite and akermanite. 
While the small deviation from ideality shown by 



Charlu et al. ( 1981 1 for the AG of mixing seems to suggest 



that the ideal solution is a good approximation to describe 
the behaviour of this phase, several studies suggest that the 
melilite phase is in fact non-ideal. Models and results from 
|Beckett| ( |1994[ ) show how melilite deviates from the ideal, 
suggesting that this behaviour has to be taken in account 
when calculating phase equilibria involving this phase. More 
recent studies by |Merlini et al.| (|2004| and |Gemmi et al.| 
( |2007[ ) further support the non-ideality of this system. 

The MELTS Supplemental Calculator was used to cal- 
culate thermodynamical properties of the melilite binary 



system using data from Charlu et al. ( 1981 1 and Walbaum 



(19731, extrapolating activities for both gehelenite and ak- 



ermanite at fixed values of temperature, pressure and mole 
fractions. We set To = 1000 K, P = 1 atm and used values 
of mole fractions from ^ Xce ^ 1 and 1 X A k ^ 0. 
Fitting polynomials to these data, we derived functions for 
the activity a and then, using the regular solution model, 
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for the activity coefficients given by 

1Ak (X Ak ,T) = (-8.98X Afc + 28.71Xl fc - 29.08Xl fe 



+8.00X 2 Ak + 2.11X Ak +0.22) 



(To/T) 



,(27) 



lGe{X G e,T) = (4.05X^ e - f 5.05X£ e + 17.75XL 

-5.97X| e - 1.24X Ge + 1.47) (T ° /T) , (28) 

for T = 1000 K. 



and fitting a polynomial, a relation between activity coeffi- 
cient and composition of spinel species at To = 1300 K was 
determined via 

lFeS P {T, X FeSp ) = (— 0.175Xp eSp + 0.373X FeSp 

+0.800) To/T , (36) 
-yM g Sp(T,X Ma s P ) = (0.035Xm 9 s p + 0.96) To/T . (37) 



3.2.7 Fassaite phase 

Activity and activity coefficient data for the fassaite phase 



were deduced by Benisek et al.| ( |2007[ |. They studied the 
thermodynamic properties of Ca-tshermak (CaAl 2 Si0 6 ) and 
diopside (CaMgSi 2 6 ), and found activity coefficient rela- 
tions 



(29) 



RTlnjo, = X 2 CaT [W m + 2{W C aT-W D% ) 
x(l-X CaT )], 

RT In 7CaT = [WcaT + 2(W Di ) ~ Wca T )X C aT] 

x(l-Xcar) 2 , (30) 
for Ca-tshermak (CaT) and diopside (Di), where W is the 



Margules parameter reported by Benisek et al. (20071 for 
both Ca-tshermak and diopside. 

The relation that links the activity coefficients and ac- 

(31) 



tivity reported by Benisek et al. ( 2007 1 is 



id 

a — a 7, 



where 



id 
a Di 
id 
a CaT 



(1- X CaT )\{2- X C aT)/2]\ 
(4XcaT)(X C aT /2)(2 - X C a T /2) . 



(32) 
(33) 



Using these relations, we can calculate activities of Ca 
Ts and diopside using a = yX together with equations ( 32 1 



and (331 for XcaT and Xd%, and equations (29 1 and (301 
for 7. 

Fitting polynomials to the calculated data, we derived 
equations for ac a T as a function of XcaT and aoi as a func- 
tion of Xr>i (where Xoi = i — Xcar)., and then for 7 accord- 
ing the regular solution model. Finally we obtained solutions 
for the activity coefficients 



(0.43X CaT - l-79X CaT + 2.32X CaT 
+0.04) To/T , (34) 
(0.17X 3 Di - 0.018X1,^ + 0.52X Dl 
+0.32) To/T , (35) 



1Co.t(T, Xcar) = 
jDi(T,X m ) = 
where T = 1500 K. 
3.2.8 Spinel phase 

The activity data for spinel species (MgAl 2 4 and FeAl 2 4 ) 
was taken from activity-composition measurements made by 
Jacob & Patill (fT998l) at T = 1300 K. They measured the 



electromotive force for different compositions of the spinel 
solution, which is linked with FeAl 2 4 activity. They de- 
duced activity-composition relationship for MgAl 2 4 using 
the Gibbs-Duhem Equation. 



3.2.9 Metal alloy phase 

Results from our ideal model simulations show that Fe 
is among the most abundant species in the solid solu- 
tion. Indeed, it represents between 90 and 94 mole-percent 
across the entire temperature range where the metal phase 
is present, with Ni representing between 6 and 10 mole- 
percent, and Si and Al less than 10 -7 mole-percent. 

The activity coefficient for Fe and Ni can be studied in 
a Fe-Ni solution, since the amounts of Si and Al are negli- 
gible and hence they will have little effect on the solution 
behaviour. Conard et al. ( 1979 1 studied the Fe and Ni be- 



haviour for the liquid-solid iron-nikel alloys for three specific 
values of temperature in the range of 1273.15-1873.15 K. 
The range of interest to us is 1273.15-1573.15 K where the 
Fe-Ni system is in solid phase. Fitting polynomials to the ex- 
perimental data, we derive the activity and then the activity 
coefficients given by 



7Fe (X Fe ,T) 
7JVi(Ajvi, T) 



lQ.27X% e + 4.854A| e 



(5.012X|, 

+1.401X Fe + 0.003) To/T , 
(-2.91X^ + 4.92X^ - 1.33Ai 
-0.27X JVl + 0.60) To/T , 



(38) 
(39) 

for To=1273.15 K. These relations are shown in Fig. [lb. 

The activity data for the silicon as solute in Fe-solution 
were taken from |Sakao fc Elliott | fl975| ). They studied the 
activity of Si in solid Fe-Si alloys in the region 0.028 ^ Xsi ^ 
0.084 and 1373.15 ^ T(K) sC 1643.15. The experimental 
results clearly show a relation between Si activity coefficient, 
T and Xsi given by 



log7s* = 1-19 



7070 



+ -6.30 



18.30 

T 



Xsi 



(40) 



Our ideal results show that Si has an amount of ~ 10~ 
mole-percent so the equation could be simplified for very 
low Si-composition as 



7si(T) = 10 



(1.19-7070/T) 



(41) 



No data for aluminium in solid Fe-Al alloys was found 
at very low Al composition, so the Al activity coefficient re- 
lationship was chosen as ideal at high temperature following 
the equation |7]l 

0.99 To/T . 



7a; (T) 
where T = 1873.15 K. 



(42) 



Using the experimental data of Jacob & Patil ( 1998 1 



3.3 Testing the regular solution model 

When we have experimental data for the same solution at 
different reference temperatures we can test the regular so- 
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0.4 0.6 
Fe Composition 

Figure 2. Test of the regular solution model for iron given by 
equation {38} using equation {H} . Here T l= 1273 K, Xj=1473 K 
and Tj.=1573 K. Polynomials of different activity coefficients are 
derived using data provided by Conard ct al. (1979). The ideal 
solution is shown as a straight line at 7 = 1. 



lution model. According to the regular solution model 

77 i -^ /7 \ (43) 



and for different experimental data at three different refer- 
ence temperatures, To = 2<,j,fc, we find 



Tj/Ti T k ITi 



(44) 



Using the set for activity coefficients by Conard et al. ( 1979) 



for iron in the metal phase for three different values of To , we 
derived activity coefficient polynomials at Tj=1473 K and 
Tk — 1573 K and tested the appropriateness of the our regular 
solution model given by equation {38} where the reference 
temperature is Ti=1273 K. Fig. [2] shows the application of 
equation (44 1. The test returns good agreement across the 



entire range of composition, with profiles given by the second 



and the third terms of equation ( 44 1 fitting the first function 
with a standard deviation of 0.03 (71573) and 0.02 (71473). 
Fig. [2] also clearly show how the behaviour of Fe deviates 
from the ideal at lower concentrations (Xf e ^ 0.5). 

We can also use experimental data to test the assump- 
tion of equation {7} . Following the results reported by Beck- 



ett 



( 1994} for gehlenite and akermanite, which show activ- 
ities moving toward the ideal with increasing temperature, 
we ran a test in which we compared results using equa- 
tions (271 and (28 1 with equation {7}. Differences between 
the results are reported in Table [4] As can be seen, the vari- 
ations in the temperature of appearance, T a , and disappear- 
ance, Td, and the maximum amount of composition, A max , 
between the results of the two models are small. Extrap- 
olations made by MELTS Supplemental Calculator (equa- 
tions (27 1 and (28l) and results given by equation {7} are in 
good agreement at the reference temperature To = 1000 K 
for the melilite phase considered here, giving us confidence 
in equation {7} when applied to species for which we do not 
have experimental data, such as NiS (equation |14[ ) and Al 
(equation |42} . 



4 SIMULATION RESULTS 

Now that we have all the thermal data, we can use our so- 
lar composition for a given pressure, P, and temperature, 
T, to solve for the equilibrium condensation sequence using 
the HSC software. As shown by |D'Alessio et al.| {1998} , the 
temperature and particularly the pressure in the disc span 
a very large range (see section [5]l). Runs were made for 
fixed values of pressure (10 -3 , 10 -4 , 10 -10 bar) over the 
temperature range 50 ^ T(K) ^ 1850. The condensation 
sequence was calculated as a function of temperature and 
results are reported using temperature-amount diagrams. 

In this section we describe the overall results of our 
equilibrium calculations for both gases and solids. Our sim- 
ulations show that the condensation temperatures depends 
on pressure, with sequences moving towards lower tempera- 
tures as the pressure decreases. This behaviour is seen more 
clearly at higer temperatures than in the cooler regions (see 
for example Fig. [3] and Fig. [4}. The temperature of appear- 
ance (disappearance) is defined as the temperature at which 
the amount of a compound raises above (falls below) 10~ r 
kmol. In our data tables, we indicate with trace those com- 
pounds which are present with an amount between 10 -9 and 
10 -7 kmol. The following discussion is made using the re- 
sults of the condensation sequence at P = 10 -3 bar unless 
otherwise specified. 



4.1 Gases 

H, H 2 and He are the most abundant gas species in the 
system and they reach a maximum amount of 8.56 x 10" 1 
kmol, 45.48 kmol and 8.89 kmol respectively in different con- 
densation regions. Looking at the entire temperature range, 
H a O, SiO, CO, N a , H 2 S, CH 4 and NH 3 are the major gas 
compounds in our results (Fig. [3}. SiO, CO and H 2 char- 
acterize the higher temperature region with CO being the 
most abundant gas. The region between 600 ^ T(K) ^ 850 
shows an important transition with CH 4 replacing CO, and 



H 2 becoming the most important gas. N 2 is present along 
almost the entire temperature range until T ~ 400 K, where 
it is replaced by NH 3 . Among the less abundant gases, sulfur 
gases (HS, SiS) are also formed. 

Moving to lower pressure, we see the entire sequence 
move to lower temperatures, including the transition be- 
tween CO and CH 4 - see Fig. [3] The maximum amount 
of each gas at equilibrium changes little with pressure. 



4.2 Solids 

The high temperature region (1400 sC T(K) sC 1850) 
is the place where the calcium-aluminium solids start to 
form, with several compounds in sequence from hibonite 
(CaAl 12 19 ) to gehlenite (Ca 2 Al 2 Si0 7 ), followed by the first 
calcium-magnesium silicates, akermanite (Ca 2 MgSi 2 7 ) 
and calcium-Tchermak (CaAl 2 SiO e ) - see Fig. [4] which 
shows the high temperature condensates for three pressures. 
The condensation temperatures of Ca(g) and Al(g) are very 
close and sensitive to the pressure. At P = 10 -3 bar, Ca(g) 
and Al(g) start to condense at T — 1738 K and the first, 
most stable solid compound formed is hibonite (CaAl 12 O ig ). 
Lowering the pressure to 10 -6 bar changes the sequence 
with the Al(g) (T = 1536 K) condensing before Ca(g) 
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Table 4. Temperature of appearance, T a , and disappearance, T^, and the maximum amount, A max , for gehlenite and akermanite calcu- 
lated by HSC using polynomials for the activity coefficients derived by MELTS Supplemental Calculator (equation | |27| l and equation ( |28| ) 
and equation {7| with T = 1000 K. 

lge = eq. g -y ak =eq. {27} 7 9e = 7 afc = 0.99 T «/ T 

(bar) (K) (K) (kmol) (K) (K) (kmol) 



Gehelenite 



io- 


-3 


1556 


1445 


9.47 x 10" 


-5 


1566 


1445 


9.43 x 10- 


-5 


10- 


-4 


1506 


1374 


1.00 x 10- 


-4 


1496 


1374 


9.52 x 10- 


-5 


io- 


-6 


1364 


1253 


9.63 x 10- 


-5 


1364 


1253 


9.63 x 10- 


-5 


io- 


-8 


1202 


1152 


9.70 x 10- 


-5 


1212 


1152 


9.70 x 10- 


-5 










Akermanite 








io- 


-3 


1556 


1445 


1.80 x 10- 


-5 


1556 


1445 


2.30 x 10- 


-5 


io- 


-4 


1496 


1374 


1.58 x 10" 


-5 


1496 


1374 


1.96 x 10- 


-5 


io- 


-6 


1364 


1273 


5.37 x 10" 


-6 


1354 


1253 


1.29 x 10- 


-5 


io- 


-8 


1202 


1152 


7.58 x 10" 


-6 


1212 


1152 


5.97 x 10- 


-6 



(T = 1506 K) and the first most stable solid compound 
formed is the corundum (A1 2 3 ). 

The central zone of our temperature range (600 ^ 
T(K) ^ 1400) is dominated by Mg-silicates, like forsterite 
(Mg 2 Si0 4 ) and enstatite (MgSiO a ), plus metal iron - see 
Fig. [5] The formation of enstatite and forsterite is chemi- 
cally related. The condensation sequence of these two com- 
pounds shows that the first peak of forsterite formed be- 
fore enstatite is pressure sensitive: lowering the pressure de- 
creases the amount of forsterite formed in its first peak until 
enstatite become the most dominant of the two (Fig.j5|. Mg- 
spincl (MgAl 2 4 ) and diopside (CaMgSi 2 6 ) also condense 
in this region. 

The cooler region (50 < T(K) < 600) shows the con- 
densation of many Fe-compounds: troilite (FeS), fayalite 
(Fe 2 Si0 4 ), ferrosilite (FeSiO s ) and Fe-spinel (FeAl 2 4 ) - 
see Fig. [6] We also see that in this temperature zone the 
condensation sequence is not very sensitive to the pressure, 
unlike the high temperature regions. From our definition of 
the temperature of appearance, which is when the amount 
of a compound reaches 10 -7 kmol, it may seem that com- 
pounds have not yet condensed, when in fact they are in 
trace amounts. For example, with the olivine and orthopy- 
roxene phases, the amount of fayalite (Fe 2 Si0 4 ) is non-zero 
in the temperature range from T — 1415 K (the tempera- 
ture of appearance of forsterite) to T = 717 K (when fayalite 
finally reaches an abundance greater than 10 -7 kmol). Sim- 
ilarly, the amount of ferrosilite (FeSi0 3 ) ranges from about 
1 x 10~ 8 to 8 x 10 -8 kmol over the temperature range from 
T = 1384 K (the condensation temperature of enstatite) 
and T = 757 K (when the amount of ferrosilite reaches 
10 -7 kmol, i.e. when it appears). 

The full list of compounds formed is given in Ap- 
pendix |B2| together with their condensation temperatures 
for different values of pressure. 



4.3 Dust distribution 

Figure [7] shows the percentage by mass of dust in the system 
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Figure 3. Evolution of the major gases as a function of temper- 
ature at (top) P = 10~ 3 bar and (bottom) P = 10 -6 bar. 



Fig. [7] shows an abrupt change between 150 and 200 K, 
which is where H 2 ice forms (see table B2 1. In a protoplane- 



tary disc, the inner edge of the region where the temperature 
falls below the condensation temperature of water is called 



for P — 10 3 and 10 6 bar. It shows that the dust mass the snow line <|Sasselov fc Lecar 20001. Our simulation at 



distribution is a function of temperature, and thus a function 
of radius, assuming that the dust temperature decreases as 
radius increases. 



P = 10~ J bar show that T = 191 K is the point where the 
process of condensation of H 2 0(g) starts, while T — 141 K 
is the temperature at which the amount of H 2 0(g) drops 
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Figure 4. Condensation sequence in the high temperature region 
as function of temperature at (top) P = 10 -3 bar, (middle) P = 
10 -4 bar and (bottom) P = 10~ 6 bar. 
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Figure 5. Condensation sequence in the middle temperature re- 
gion as function of temperature at (top) P = 10 — 3 bar, (middle) 
P = 1Q- 4 bar and (bottom) P = 10~ 6 bar. 



below 10 -7 kmol. At P = 10 -6 bar, the temperature in 
which the process of vapor condensation starts is T = 151 K 
while T — 110 K is the temperature at which the amount of 
H 2 0(g) drops below 10~ 7 kmol. 



5 REGULAR VERSUS IDEAL SOLUTION 
MODELS 

In order to understand the differences between the regu- 
lar and the ideal solution models, we need to return to the 



definition of activity, a, and activity coefficient, 7. Activity- 
composition diagrams derived from experimental data pro- 
vide information about the departure from the ideal be- 
haviour of compounds according to their composition at 
the reference temperature (see Fig. [T]). Activity coefficient- 
composition diagrams, derived from equations in sect ion [372] 
also provide information about the departure from the ideal 
with varying temperature across the entire composition 
range of a selected compound. Figures [8] and [9] show ac- 
tivity coefficient polynomials for forsterite and enstatite at 
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Figure 6. Condensation sequence in the low temperature region 
as function of temperature at (top) P = 10 — 3 bar, (middle) P = 
10 -4 bar and (bottom) P = 10~ 6 bar. 



different values of temperature. The 7 polynomials curves, 



from equations (18 1 and (19), move closer the ideal (7 = 1) 
with increasing temperature. This behaviour is seen for ev- 
ery compound in the system. Furthermore, every 7 curve 
for different temperatures return values close to one when 
X — > 1 for all temperatures. 

Equation Q for the Gibbs free energy of compounds at 
a given temperature could also be written as 

Ad = Gi-G° = RTlaji+RTlnXi = AGf+AGf, (45) 

where RTlaXi is the partial molal Gibbs free energy of 
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Figure 9. Activity coefficient of enstatite (MgSiO a ) across the 
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equation JTsJ . The 7=1 line represents the ideal. 
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mixing for an ideal solution, and RTlnji is the departure 
from ideal behaviour of component i in the mixture. We 
define the former as AG; d and the latter AGf x , where the 
superscripts id means ideal and ex is excess. We can study 
the non-ideal contribution in the co mput ation of AGi by 
plotting AGT + AG] d from equation ( 45 1 compared to the 
ideal AG\ d at fixed values of temperature. This is shown in 
Fig. [To] for forsterite and enstatite, and we find that AG1 X 
could be high enough to change the value of the AGi. 

We compare the condensation sequences calculated us- 
ing the ideal solution and the regular solution to show how 
the activity coefficients can affect the results. Fig. |ll| shows 
the condensation sequences with both the ideal solution and 
regular solution models at a pressure of 10" 3 bar. Similar 



to define the list of possible compounds and the activity- 



plots (not shown) were obtained for P = 10 . In tables Bl 
and |B2| we also report the temperature of appearance and 
disappearance for the ideal and regular solution models. 

We see some differences in the values of the temper- 
atures of appearance and disappearance. The first peak of 
forsterite at 1394 K is more abundant in the regular solution 
than for the ideal case at P = 10 -3 bar. In the low tempera- 
ture region, the regular solution model returns a larger range 
of temperatures in which fayalite (Fe 2 Si0 4 ) is stable. In the 
ideal solution, this range is smaller and there is more Fe 
available for the the formation of FeS, Fe metal, ferrosilite 
(FeSi0 3 ) and Fe-spinel for T > 400 K, and at very low tem- 
perature magnetite (Fe 3 4 ) replaces fayal ite (this also hap 
~ 6 bar) 



11 



Furthermore, 



pen at the pressure of 10 bar) - see Fig. 
there is no evidence of formation of magnetite (Fe 3 4 ) in 
the regular solution at the temperature and pressure range 
used in these calculations (see Fig. Rjl. 

The ideal solution model also shifts the turning point of 



enstatite and forsterite to lower temperatures (see Fig. Ill, 
which also explains the lack of fayalite (see section [6. 1.1[ ) in 
these temperature regions. 

Note that the errors in temperatures are 5 K due to the 
step size used in the HSC calculations. 



6 DISCUSSION 

In this section we compare the results our regular solution 
model with previous thermodynamic models and observa- 
tional data. 



6.1 Comparison with models 

In order to check the results of our regular solution model, 
we compare our simulation results to those of earlier studies 



by Pasek et al. ( 2005 1, Yoneda & Grossman ( 1995 I and|Gail 



(19981. The main difference between our model and theirs 



is the regular solution model we described in section [3] 



Pasek et al. ( 2005 \ carried out their calculation using 



the HSC software package, using 150 gas species and 100 
solid species from the 20 most abundant solar elements, but 
they assumed their solid solutions as ideal. Their calcula- 
tions assumed a temperature range of 400 ^ T(K) ^ 2000. 

|Yoneda fe Grossman| ( |1995[ ) created a larger system 
with solid, liquid, and gas states. They used 120 gas species, 
5 liquid species, and 110 solid species from the 20 most abun- 
dant solar elements. They studied the behaviour of the liquid 
solution and each solid phase using data from the literature 



composition model from Berman (1983) for non-ideal liquid 



solution. They calculated equilibrium for temperatures from 
900 K to 1800 K, and with a larger range of pressures (from 
10~ 6 to ~ 1 bar). Their results show that the liquid phase 
appears only at high pressures (P ^ 0.1 bar), confirming 
the solid-gas assumption in the pressure range used in our 
model. 



Gail ( 1998 1 studied the chemical equilibrium between 



the gas phase and the solid species, and determined rela- 
tionships between the partial pressure of the gas involved in 
the creation of the solid species and the equilibrium constant 
of the reaction. 

Our calculations are in good general agreement with 



Pasek et al. (20051, Yoneda & Grossman (19951 and Gail 



( 1998 1 in terms of the condensation sequence and the tem- 



perature of appearance and disappearance. However, there 
are some important differences. Table |C1| shows our con- 
densation temperatures compared with these three models. 
Small differences between the temperature of appearance 
and disappearance can be seen, except for the plagioclase 
phase for which large discrepancies are seen, and the order 
of appearance of forsterite and enstatite in the low pressure 
regime. However, our regular solution model and our ideal 
solution model (see tables Bl and B2| are consistent. We 



will discuss forsterite and enstatite and the discrepancy in 
the plagioclase in more details in the next section. 

Overall, the condensation sequence is the same for the 
12 most abundant species in our model and those in |Pasek] 
et al. ( 2005 1 and Yoneda & Grossmanl (119951) except for the 



corundum (A1 2 3 ). Indeed, the previous models show that 
the corundum is created just before the hibonite species, 
whereas in our model, hibonite is the most refractory species 
of the system. This difference between the models comes 
from the closeness of the melting point between Ca and Al. 



In the model of Yoneda & Grossman (19951, Ca begins to 



condense in hibonite at a lower temperature than the Al 
melting point, and thus A1 2 3 can be created between these 
two condensation temperatures, whereas our system does 
not show these difference of melting point temperatures un- 
til the pressure drops to P — 10 -6 bar. Gail ( 1998 \ reports 



corundum as the first solid formed, followed by melilite. The 
absence of hibonite (CaAl 12 O ig ) in his model does not allow 
us to compare the difference of our corundum and hibonite 
condensation sequence, however the condensation tempera- 
ture of corundum at P = 10 -6 bar show similar values. 



6.1.1 Enstatite and forsterite 



The behaviour of forsterite and enstatite see n by |Yoneda 
& Grossmanl (l995l), [Gaul (l998|), |Ganl |2004l, iPasek et al 



( 2005 1 and this work all show differences from one another. 



Yoneda & Grossman ( 1995 1 report in their Fig. 2 an ini- 



tial peak of forsterite at high temperatures. The abundance 
of forsterite then drops and enstatite replaces it as the most 
stable compound at T = 1366 K. There is then a plateau of 
constant enstatite/forsterite ratio until T = 1188 K where 
the amount of forsterite increase while enstatite remains 
constant. Yoneda & Grossman ( 1995 1 have T — 900 K as 



the lower limit of their calculations and so it is not possible 
to compare results in the lower temperature regions. 
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Figure 11. Condensation sequences in the low temperature region (left), the middle temperature region (middle) and the high temper- 
ature region (right) for the regular (top) and the ideal (bottom) solution models. Pressure is 10 — 3 bar. 



Gail (1998) shows the condensation sequence of en- 



Pasek et al. ( 2005 I show the equilibrium condensation 



statite and forsterite in terms of radial variation of the frac- 
tion of silicon condensed into these compound. The first 
peak of forsterite is followed at T — 1319 K by an in- 
crease of enstatite that becomes the most stable compound 
- see their Fig. 21. Then there is a plateau of constant 
enstatite/forsterite ratio and there is no further change 
in either enstatite or forsterite (out to R = 2 AU and 
T sC 300 K). 



In Gail ( 2004 1 an initial peak of forsterite precedes the 
formation of enstatite. Then enstatite overtakes forsterite 
at 0.55 AU (T ~ 1350 K) - see their Figs. 10 and 11 - and 
then rather than a plateau of constant enstatite/forsterite 
ratio, both the enstatite and forsterite abundances drop with 
decreasing temperature. 



sequence in terms of atom fraction of cationic elements. Dif- 
ferences from previous work can been seen in their Fig. 2, 
in which the amount of enstatite continues to increase once 
fayalite forms at T ~ 450 K. Our condensation sequences 
shows a complex behaviour, demonstrating that forsterite 
and enstatite follow a condensation sequence of forsterite- 
enstatite-forsterite, going from the warmer region to the 
cooler outer region of the disc (see Fig. [5J . Forsterite first 
condenses at T = 1415 K and is then replaced by enstatite 
(with has a temperature of appearance of T = 1384 K ). As 
the temperature decreases (1000 ^ T(K) ^ 1300), a plateau 
of constant enstatite/forsterite ratio can be seen. At T ~ 
1000 K there is a sudden drop of enstatite and an increase 
in the abundance of forsterite. A second plateau of constant 
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enstatite/forsterite ratio is seen through 600 ^ T(K) ^ 900. 
At T ~ 600 K forsterite increases its abundance, overtak- 
ing enstatite at T ~ 550 K, while the enstatite abundance 
continues to drop with decreasing temperature. The tem- 
perature values reported here refer to a pressure of P=10~ 3 
bars. 

The differences in the condensation sequence reported 
by these four references and our work could be due to several 
factors, such as different initial compositions, different lists 
of possible compounds considered in the final system, differ- 
ent solution model and different thermodynamical databases 
used in the calculations. To resolve this issue, which is of 
great importance to astrophysics and in understanding ob- 
servational evidence from both protoplanetary disks and 
meteorite data, a controlled laboratory experiment using a 
complex system over a range of pressures and temperatures 
is required. 

While we cannot explain the differences in the forsterite 
and enstatite behaviour seen in these four references, we can 
explain the forsterite and enstatite behaviour seen in our 
results. Here we will discuss (i) the formation of enstatite 
after the first peak of forsterite, (ii) the drop of enstatite and 
the increase of forsterite at T ~ f 000 K, and (iii) the turning 
point between enstatite and forsterite around T ~ 600 K. 



Schegerer et al. 



1 2006 ) showed that the reaction from 
forsterite to enstatite is possible if quartz (Si0 2 ) is present 
according to the following reaction given by|Rictmcije r et al.| 
(19861: 



Mg 2 Si0 4 + SiQ 2 -» 2MgSi0 3 



(46) 



In our simulation, there is no Si0 2 present as an isolated ox- 
ide. This is because simulations made using the Gibbs free 
energy minimization do not provide information about the 
rate of the reaction and the steps between the initial and the 
final conditions. In our equilibrium conditions, silica is al- 
ways involved in the formation of more complex compounds. 
The amount of SiO(g) in our simulation drops by four orders 
of magnitude in the range of temperature in which the two 
silicates form (1152 sC T(K) sC 1450). 

The formation of enstatite from forsterite at higher tem- 
perature could be also chemically explained with reactions 
involving Mg(g) with the condensation of SiO(g): 

Mg 2 Si0 4 + O(g) + SiO(g) ->■ 2MgSi0 3 , (47) 

but this reaction is related to the availability of oxygen and 



to the kinetic processes that drive its rate (see Pegourie & 
Papoular 1985 and reference therein). Given the low free 
oxygen abundance assumed to exist in protoplanetary disks, 



a more likely reaction is the following ( Imae|1993 1: 

Mg 2 Si0 4 + H 2 0(g) + SiO(g) -> 2MgSi0 3 + H 2 (g). (48) 

This reaction is also supported by the decrease of H 2 0(g) 
and SiO(g) we see in Fig. [3] in the high temperature region. 

In our model we find that when T ~ 1000 K, the amount 
of enstatite drops dramatically and forsterite starts to form 
again. 

This sudden drop of enstatite can be explained with the 
following three reactions that occur simultaneously in this 
temperature region: 



2MgSi0 3 + CaMgSi 2 6 + MgAl 2 4 
CaAl 2 Si 2 8 + 2Mg 2 Si0 4 , 
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Figure 12. Condensation sequence from T = 500 K to T = 
1200 K for P = 10 -3 bar for compounds involved in the reac- 
tions | |49| l [top] , ( |50| [middle] and l |51| [bottom] using the regular 
solution model. The amount of anorthite, Mg-spinel, diopside, al- 
bite, Ca-Ts and Na(g) are multiplied by a factor of 10 for clarity 
in this plot. 



from Marfunin (19981, 



(49) 



13MgSi0 3 + MgAl 2 4 + 2Na( 5 ) + H 2 0(ff) -» 
7Mg 2 Si0 4 + 2NaAlSi 3 8 + H 2 (.g), 

from |Lewis"| ( |2004[ ), and 

Mg 2 Si0 4 + CaAl 2 Si 2 O g 



(50) 



(51) 
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Figure 13. Fraction of Si-bearing compounds in the 500— 
1200 K temperature range, where reactions ( |49[ ), ( |50[ |, 
( |51[ ) and | |54| | take place. 



from Marfunin (19981 (see Fig. 121. The forsterite:enstatite 
mole ratio produced by these three reactions is 10/17 = 0.59. 
We can compare this forsterite:enstatite ratio with the data 
returned by our calculation. We define AX en and AXf a by 



AA en 
AX fo 



= IX 
= \x 



So 



(52) 
(53) 



which represent the difference between the amount of com- 
pound at T = 1000 K and T = 950 K for enstatite, 
AX en , and forsterite, AXf a respectively. Our calculation 
returns AX en =0.0010477 and AX fo =0.000593. Thus 
AX f / AX en = 0.56, in good agreement with the mole ratio 
involved in these three reactions above. 



According our simulations, reactions ( 49 1, ( 50 \ and ( 51 \ 



become most favourable when the temperature reaches T ~ 
1000 K. These reactions also explain the difference between 
our temperature of appearance of plagioclase (T = 980 K) 
and those derived by previous works. 

In the mid-temperature range 650 ^ T(K) ^ 900 the 
enstatite-forsterite ratio is constant, with T ~ 550 K be- 
ing the turning point where forsterite increases its amount 
over enstatite (see Figs. [5] and |12[ ). Thermodynamic pre- 
dictions by Fegley (2000) involving Fe show that, with de- 
creasing temperature, the following reaction becomes more 
favourable: 



2Fe + 2MgSi0 3 + 2H 2 0(g) -» 
Mg 2 Si0 4 + Fe 2 Si0 4 + 2H 2 (ff). 



(54) 



This reaction explains the conversion of enstatite to 
forsterite at T ~ 600 K (the turning point) and also the 
formation of fayalite (Fe 2 Si0 4 ) at lower temperature (see 
Figs. |6| a nd JT2| . 

Fig. |13| show the Si-bearing in the 500-1200 K tempera- 
ture region in which reactions ( ]50[ ), ( |51[ ) and | ]54| take 
place. 

In terms of maximum abundance, we find that enstatite 
is the dominant compound among the crystalline silicates, 
as is demonstrated in previous works. However, our model 
suggests this is only the case in a well determined range of 
temperatures going from T ~ 600 KtoT~ 1400 K. 

As briefly introduced in Section |4.2| according to our 
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Figure 14. Difference between the temperatures of appearance 
of forsterite and enstatite (AT) at different pressure values as 
reported by |Yoneda fe Grossman| ( |1995^ (blue diamonds) and 
this work (red dots). (See Table|5(for data values from this work) . 
The AT = line represents the point in which T a (forsterite) = 
T a (enstatite). 



calculations made using HSC with the regular solution 
model, the first peak of forsterite is pressure sensitive (see 
Fig. pk. We also see that the difference between the tem- 
peratures of appearance of forsterite and enstatite becomes 
smaller as the pressure decreases. When the pressure drops 
to P = 10~ 6 bars, HSC returns T a = 1212 K as the tem- 
perature of appearance for forsterite and T a = 1223 K for 
enstatite. This in contrast with previous results reported 
by |Yoneda fe Grossman| ( [T995] , [Gail] ( p98"] ), [Gall] pOQ4| 
andlPasek et al. ( 2005 1. Table |5 reports the temperatures of 



appearance of forsterite and enstatite at different values of 
pressure (from P = 10 -3 bars to P = 10~ 6 bars) in our cal- 
culations. The temperature difference, AT = T a /or - T a ens , 
is also shown. 

The decreasing difference between the temperatures of 
appearance of forsterite and enstatite as the pressure de- 
creases was also noticed by Grossman (19721. Indeed, the 
data reported by |Yoneda fe Gr ossman ( 1995 ) shows a similar 
behaviour when the press ure decreases from P = 10~ 3 atm 

14 we plot this temperature dit- 



to P = 10~ b atm. In Fig. _ 
ference, AT, with pressure as found by Yoncda & Gro ssman] 
( 1995 1 and in our work. Similar results were also seen by 



Gail 



Gail 



( 1998 1 and|Gail| (|2004||. Yoneda fe Grossman] (|1995[ ), 
( T998| , |Gail| ( |1998| |2004| > and |Pasek et al.| ( |2005| ) do 

not show the swap in the order of of appearance of forsterite 
and enstatite at P = 10 -6 bar which we find with our reg- 
ular solution model. We can only assume that it is due to 
either the solution models used, or more likely, the thermo- 
dynamics data used or indeed a combination of both. 



6.2 Comparison with observations 

6.2.1 Discs 



|Pollack et al.| ( |1994[ ) studied the mass fraction of the eight 
most abundant elements (H, C, O, N, Mg, Si, Fe, and S) in 
molecular clouds, the interstellar medium, and solar system 
bodies. They also deduced a model for the composition of 
each compound for the outer parts of protoplanetary discs 
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Table 5. The temperature of appearance, T a , of forsterite and 
enstatite at different values of pressure, P , and the difference in 
their temperatures of appearance, AT. 



P (bar) 


Ta (K) 


T a (K) 


AT (K) 




Forsterite 


Enstatite 
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1344 


1334 
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icr 5 
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1278 


7 


icr 6 


1212 


1223 


-11 



where T ^ 700 K. Table[6]shows the distribution of elements 
in a typical disc, considering three kinds of matter: ices, 
refractory grains (which condense at high temperature), and 
gas. The main refractory elements observed by |Pollack et al.| 
( 19941 are the same as in our model: olivine, orthopyroxene, 



ment and support our condensation results. Our equilibrium 
simulation shows that forsterite becomes more stable than 
enstatite in the low temperature region (150 ^ T(K) ^ 500), 
with enstatite more abundant when T 500 K. 



6.2.2 Meteorites 

Additional information can be inferred from the study of 
chondrites. Chondrites are thought to be the oldest known 
components of the solar system and are assumed to have 
been created during the protoplanetary disc phase. Thus, 
they provide important evidence of the conditions present in 
the early stages of the solar system's formation. Chondrites 
are composed of three major parts: Calcium and Aluminium 
Inclusions 



troilite, and Fe species. Indeed these species are present for 
a large scale of temperature and for a high mol amount. 



The results derived by Pollack et al. (19941, based on 



astronomical observations and modeling, show that most of 
the silicates are present in two main phases, olivine and 
orthopyroxene, with the former dominating. In the mid- 
dle temperature region of our condensation sequence (see 
Fig. pjl, moving toward lower temperatures we see that 
forsterite, the main olivine species in our system, replaces 



(CAIs), chondrules and matrix (Scott 20071. 
CAIs are the oldest objects in the chondrites and they 
are created at highest temperatures. CAIs are composed 
of different phases, including spinel (MgAl 2 4 ,FeAl 2 4 ), 
melilite (Ca 2 Al 2 Si0 7 ,Ca 2 MgSi 2 7 ) and compounds like 
perovskite (CaTiO s ), hibonite (CaAl 12 19 ), calcic py- 
roxene (CaMgSi 2 O s ,CaAl 2 SiOg ,CaTi Al 2 O e ,CaTiAlSiO e ) , 
anorthite (CaAl 2 Si 2 8 ), and forsterite (Mg 2 Si0 4 ), and also 
include less common compounds like grossite (CaAl 4 7 ), 
Ni-Fe metal, corundum (A1 2 3 ), rhonite (Ca 4 (Mg,Al,Ti)i2- 
(Si,Al)i 2 O 40 ) and CaAl 2 4 ( |MacPherson|2006[ >. 

All the most important elements are present in our sys- 
tem at high temperatures (T > 1300 K at P = 10~ 3 bar). 
This first comparison provides a validation of the condensa- 
tion process. 

Chondrules are composed of less refractory species and 
their structure is very diverse. However most chondrules con- 
tain the same main species of olivine (Mg 2 Si0 4 , Fe 2 Si0 4 ), 
pyroxene (MgSiO s , FeSi0 3 ), metallic Fe-Ni and Fe-Ni sul- 
perature for troilite derived with our regular solution model fides (FeS, Ni 3 S 2 , NiS). The results of our regular solution 

model support these observations, as these species are the 
main Fe, Mg, Ni compounds in our model creating in the 
middle and low temperatures regions (T ^ 1400 K with a 
pressure of 10~ 3 bar). 

The high temperature region between 1850 K and 
1400 K provides information about the first solids that con- 
densed in protoplanetary discs. It is the region in which Ca- 
Al compounds start to condense into a variety of species. Ac- 
cording our model, the condensation temperatures for Al(g) 
and Ca(g) are sensitive to the pressure (as shown in Fig.[4|. 
At higher values of pressure, when Ca(g) condenses at the 
same temperature as Al(g), no corundum (A1 2 3 ) forms as 
an isolated compound. This is because, at equilibrium the 
corundum is directly involved in the formation of Ca-Al 
compounds like hibonite according to the reaction 



enstatite as the most stable compound. Pollack et al. ( 1994 1 
indicate troilite (FeS) as the most abundant Fe-compound 
(T $J 680 K), also stating that its temperature of conden- 
sation is almost insensitive to the pressure. Furthermore, as 
the temperature increases, troilite reacts with molecular hy- 
drogen forming H 2 S(g) and Fe(s). Our results are in good 



agreement with Pollack et al. ( 1994 1. The condensation tern- 



is T = 687 K (see Tab. B2 I over a wide pressure range and 
H 2 S(g) is the most abundant sulfur-gas when the tempera- 
ture raises above T = 650 K - see Fig. [3] 

Silicates are observed almost everywhere that dust can 
survive: in the interstellar medium (ISM), in protoplane- 
tary discs, in the Earth's mantle, in comets and meteorites 
in different structures and shapes. Infrared spectral analy- 
sis shows that the ISM is characterized by amorphous sil- 
icates, while in protoplanetary discs there is evidence of 
crystalline silicates, clear proof of processing and growth 
of grains ( jOlofsson et al.||2009"| |Meeus et al.||2009[ ). Among 



them, forsterite (Mg 2 Si0 4 ) and enstatite (MgSi0 3 ) are the 
best tracers of the condensation processes involved in the 
protoplanetary dies because of their related chemistry (sec- 
tion [67TTTj) . iBouwrnatTetjiL] |2008]) used infrared (5-35 um) 
spectra to study silicates in seven protoplanetary discs 
around pre-main sequence stars with Spitzer. This wave- 
length range probes both warm (A ~ 10 um, T ~ 500 K) 
and cool (A ~ 30 um, T ~ 120 K) region of the disc. They 
found variations in the abundance of the two main silicates 
with enstatite dominating the warmer region of the discs 
and forsterite the colder region. Meeus et al.| ( |2009 1 also 
studied infrared spectra (7.5-35 urn) of 12 T Tauri systems 
and report the forsterite-enstatite ratio {Rtfj e {) for seven 
objects. The warmer region (7.5-17 urn) is characterized 
by R(f/e) ~ 0.3 while the cooler region (17-37 um) has 
0.9 < R(f/ e ) ^ 1 -2. the observation and modeling by|Bouw- 



CaO(g) + 6A1 2 3 -> CaAl 12 19 . 



(55) 



Lowering the pressure to P = 10 provides a small temper- 
ature window in which there is no CaO(g) available for re- 
action and hence Al 2 O s is the first solid compound. As soon 
as Ca(g) condenses, all the corundum is then involved in 
the formation of hibonite. Observational evidence of isolated 



man et al. (2008) and Mccus ct al. (2009) are in good agree- 



corundum has been reported by Nakamura et al. (20071. 
They found forty-three corundum grains (1-11 um in size) 
in the matrix of the carbonaceous chondrite Acfer094 by 
using cathodoluminescence imaging. Their petrographic ob- 
servations indicate that, in the framework of the homoge- 
neous condensation theory, the pressure had to be lower 
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Table 6. Fractional abundance of elements between gas, ice, and refractories in protoplanetary discs jPollack et al.|1994| 



Elements Gas Refractory Ice 



Si 1 

(Fe 2 Si0 4 , Mg 2 Si0 4 , FeSi0 3 , MgSi0 3 ) 

Mg 1 

(Mg 2 Si0 4 , MgSiO a ) 

Fe 1 

(Fe 2 Si0 4 , FeSi0 3 , FeS, Fe) 

S 0.20 0.75 0.05 

(H 2 S) (FeS) (H 2 S, SO, S0 2 ) 

O 0.13 0.28 0.59 

(CO, O, H 2 0) (Silicates, CHON) (H 2 0, CO, Volatile Organics) 

C 0.28 0.55 " 0.17 

(CO, C) (CHON) (CO, C02, CH 4 ) 

N 0.79 0.20 
(N 2 , N, NH 3 ) CHON 



than 3 x 10 -5 bar. They calculate the homogeneous con- 
densation temperature of corundum as a function of total 
pressure and cooling time, showing that at pressures higher 
than few times 10~ 5 bar, the homogeneous condensation 
temperature of corundum is lower than the equilibrium con- 
densation temperature of hibonite. This means that at these 
conditions hibonite is the first compound formed. The ev- 
idence of corundum grains in Acfer94 constrains the total 
pressure in the environment in which they formed to be lower 
than 10 bar. Thus corundum, as an isolated compound, 
could be tracer of low pressure environment. 



7 CONCLUSIONS 

Determining the condensation sequence and the physical 
chemistry that drives it is an important step in understand- 
ing the first stages of the planet formation process. Previ- 
ous work made use of the ideal solution model (a = X) 
to describe the thermodynamic systems. However, labora- 
tory experiments show that the behaviour of most species is 
non-ideal. In this work we investigated the thermodynamic 
properties and the activities of several solid compounds in 
a protoplanetary discs environment specifically forming on 
the regular solution model. Many compounds in our sys- 
tem deviate from the ideal across the entire range of com- 
position, while others may be relatively close to ideal for 
most of the composition range. Decreasing the temperature 
causes further deviations from the ideal and, at fixed values 
of temperature, the composition of compounds can also vary. 
As a consequence, the activity coefficients affect the results 
of the computation of the Gibbs free energy minimization. 
From the experimental data describing the behaviour of the 
activity of several compounds in solution, and the regular 
solution model we derived, it is clear that a non-ideal solu- 
tion must be taken in account outside the very high temper- 
ature regions (where most compounds are close the ideal) 
and when the compounds' composition in the phase is not 
dominant (i.e. when X is low). We introduced activity coef- 
ficient relations for nine phases and twenty-one compounds. 
These relations can easily be used in Gibbs free energy min- 
imization tools requiring polynomials describing r yk{Xk,T). 
Our regular solution model shows that the solids which 



condense in protoplanetary discs are mainly composed of 
silicates for a large range of temperatures and of water ice 
for the cooler regions of the disc. Calcium-aluminum com- 
pounds and iron compounds characterise the high tempera- 
ture region and the low temperature region respectively. 

Our results also demonstrate that pressure plays an im- 
portant role determining the temperature of appearance and 
disappearance of compounds, their amounts (e.g. the peak of 
forsterite in the middle temperature region) and formation 
(e.g. corundum in the high temperature region). 

This work provides a starting point for the study of the 
non-ideal behaviour and the resulting composition of grains 
in protoplanetary discs. In future work, we will introduce ID 
and 2D disc models to study in more detail the condensation 
processes and the resulting condensation sequence in well- 
defined locations in the disc. 
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Table Al. Full list of species used in the models. 
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Table Bl. Temperature of appearance and disappearance for solid compounds in the system at different pressures for the ideal solution 
model. The temperature of appearance (disappearance) is defined as the temperature in which the amount of compound raises above 
(falls below) 10~ 7 kmol and we indicate with trace, compounds that are present in the results with an amount of 10~ 9 ^ kmol ^ 10~ 7 
and with <50, compounds present with an amount not less than 10 -7 kmol at the lower limit of our range of temperature (50 K). 
Compounds that do not appear in our calculations are indicated with -. 
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Table B2. Temperature of appearance and disappearance for solid compounds in the system at different pressures for the regular solution 
model. The temperature of appearance (disappearance) is defined as the temperature in which the amount of compound raises above 
(falls below) 10~ 7 kmol and we indicate with trace, compounds that are present in the results with an amount of 10 -9 ^ kmol ^ 10~ 7 
and with <50, compounds present with an amount not less than 10 -7 kmol at the lower limit of our range of temperature (50 K). 
Compounds that do not appear in our calculations are indicated with -. 
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